The many dimensions of combination therapy: How to combine antibiotics to limit resistance evolution

Abstract In combination therapy, bacteria are challenged with two or more antibiotics simultaneously. Ideally, separate mutations are required to adapt to each of them, which is a priori expected to hinder the evolution of full resistance. Yet, the success of this strategy ultimately depends on how well the combination controls the growth of bacteria with and without resistance mutations. To design a combination treatment, we need to choose drugs and their doses and decide how many drugs get mixed. Which combinations are good? To answer this question, we set up a stochastic pharmacodynamic model and determine the probability to successfully eradicate a bacterial population. We consider bacteriostatic and two types of bactericidal drugs—those that kill independent of replication and those that kill during replication. To establish results for a null model, we consider non‐interacting drugs and implement the two most common models for drug independence—Loewe additivity and Bliss independence. Our results show that combination therapy is almost always better in limiting the evolution of resistance than administering just one drug, even though we keep the total drug dose constant for a ‘fair’ comparison. Yet, exceptions exist for drugs with steep dose–response curves. Combining a bacteriostatic and a bactericidal drug which can kill non‐replicating cells is particularly beneficial. Our results suggest that a 50:50 drug ratio—even if not always optimal—is usually a good and safe choice. Applying three or four drugs is beneficial for treatment of strains with large mutation rates but adding more drugs otherwise only provides a marginal benefit or even a disadvantage. By systematically addressing key elements of treatment design, our study provides a basis for future models which take further factors into account. It also highlights conceptual challenges with translating the traditional concepts of drug independence to the single‐cell level.


| INTRODUC TI ON
Antibiotic resistance poses a major challenge for patient treatment worldwide.The rapid evolution and spread of resistance results in a loss of treatment options, increasing the demand for new antibiotics.
Designing treatment strategies that limit the evolution of resistance can potentially extend the lifetime of the antibiotics currently in use and therefore secure the availability of effective treatments for the future (Tyers & Wright, 2019).One possibility to reduce the risk of resistance evolution during treatment is to increase the genetic barrier to resistance, that is, to increase the number of mutations (or resistance genes) the bacteria need to escape the treatment.This can especially be achieved by using more than one drug in combination.Ideally (i.e. in the absence of cross-resistance), bacteria then need to acquire multiple mutations-one for each drug-to become fully resistant as opposed to a single mutation, which would be sufficient for resistance to monotherapy.The most well-known example of the efficiency of combination treatment of bacterial infections is the treatment of Mycobacterium tuberculosis infections, for which drug combinations have been the standard treatment for more than 70 years due to their ability to reduce resistance evolution (Fox et al., 1999;Medical Research Council, 1950).As drug resistance is becoming a problem for the treatment of more and more different bacterial infections, many researchers in evolutionary medicine highly recommend using combination therapy more widely (Andersson et al., 2020;Bonhoeffer et al., 1997;Merker et al., 2020;Palmer & Kishony, 2013;Roemhild & Schulenburg, 2019;Woods & Read, 2023).
In order to optimally design combination treatments, we need to gain a thorough understanding of how the choice of drugs and their precise administration affect the evolution of resistance.To prevent the accumulation of resistance mutations up to full resistance, a good combination treatment must limit the appearance of new mutations and the spread of partially resistant genotypes.
The efficiency of combinations in limiting resistance evolution depends on drug characteristics such as the bacteriostatic versus bactericidal activity (in the following referred to as modes of action), their dose-response curves, interactions between the drugs and collateral effects and complex physiological responses such as hysteresis.Interestingly, the effects of drug-drug interactions and collateral effects-two rather complex factors-on resistance evolution are more frequently explored (Barbosa et al., 2018;Hegreness et al., 2008;Munck et al., 2014;Pena-Miller et al., 2013) than those of the modes of action and the shape of the dose-response curves, which seem more 'fundamental'.For comprehensive reviews on the former, see, for example, Bollenbach (2015), Baym et al. (2016) and Roemhild et al. (2022).Several studies explored the effect of the drugs' modes of action on the type of drug-drug interactions-synergism or antagonism- (Chandrasekaran et al., 2016;Lázár et al., 2022;Ocampo et al., 2014).Coates et al. (2018) further showed that the efficiency of a combination in clearing a susceptible population depends on the modes of action of the combined drugs.Yet, to our knowledge, a direct assessment of how the modes of action of drugs in combination influence resistance evolution is missing.Besides the mode of action, the shape of the dose-response curve is another factor to consider when combining drugs.A study by Russ and Kishony (2018) showed, for example, that drugs with steep dose-response curves need to be applied at a higher total drug concentration to inhibit the growth of a bacterial population to the same extent as monotherapy.
For optimal combination therapy, we do not only need to know which drugs to apply but also how.With multiple drugs in combination, applying the different drugs, each at the same concentration as in monotherapy, may not be necessary to clear the wild-type infection.In fact, lowering the doses might even be required to avoid toxicity.Increasing the number of drugs combined while keeping the total drug dose constant can reduce the drugs' effect on the growth of the bacterial population (Lázár et al., 2022;Russ & Kishony, 2018).
How such a reduction in individual doses would impair the efficiency of combination therapy in limiting resistance evolution is unclear, especially with a larger number of drugs.While the above studies split the total drug dose equally among the drugs, unequal drug ratios could be beneficial, especially when the drugs vary in their characteristics.
In this study, we set up a stochastic pharmacodynamic model to test the efficiencies of various combination treatments when resistance evolution relies on mutations in the chromosome.We compare treatments based on their probability of successfully eradicating a bacterial population, which we derive from branching process theory.
Besides calculating the probability of treatment success, we also look at the two factors that determine the risk of treatment failure due to resistance evolution-the expected number of resistant mutants and their chance to escape stochastic loss while rare.This allows us to understand why one treatment works better than another.We focus on the effects of the modes of action, the dose-response relationships, the number of drugs and their doses and exclude drug-drug interactions and collateral effects.Calculating the combined effect of drugs in the absence of drug-drug interactions requires choosing a null model for drug additivity.In our analysis, we employ the two most common models, Bliss independence (Bliss, 1939) and Loewe additivity (Loewe & Muischnek, 1926).Bliss (1939) calls two drugs independent if the fraction of cells surviving combination treatment is the product of the single-drug survival probabilities.Loewe and Muischnek (1926) consider two drugs independent when a fraction of one drug can be exchanged by a fraction of the other without a change in drug effect.
Our systematic analysis provides an overview over the influence of a range of treatment choices and allows us to formulate some rules of thumb for the design of promising treatments.

| General model
We consider a bacterial population of initial size N 0 , undergoing treatment with n drugs.Bacteria can be fully susceptible wild-type cells (denoted by W) or resistant to any subset I of drugs (denoted by M I with I ⊆ {1, 2, … , n}).A bacterium acquires a resistance mutation to drug i during replication with probability u i .We assume that each mutation confers resistance to exactly one drug, there is hence no cross-resistance.Multiple mutations can appear simultaneously.The initial population is either fully susceptible, or resistant cells pre-exist at low frequency prior to treatment.Whenever we include standing genetic variation, we assume that the population is at mutation-selection balance (see Appendix A).We ignore resource competition and other interactions between bacteria and assume that cells replicate and die independent of each other.The population dynamics are thus described by a multi-type branching process.
In the absence of treatment, wild-type bacteria replicate at an intrinsic rate W 0 and die at an intrinsic rate 0 , resulting in a net per capita growth rate W (0) = W 0 − 0 .We assume that resistance to drug i entails a cost, reducing the replication rate by a factor 1 − i , hence 0 for a single mutant with resistance to drug i .The cost is multiplicative for types that are resistant to multiple drugs.
During treatment with drugs at concentration c 1 , c 2 , ⋯ , c n , the net per capita growth rate X c 1 , c 2 , ⋯ , c n of type X (with ing the effect of the antibiotic on the growth rate: In the following sections, we provide derivations of the drug effect E X c 1 , c 2 , ⋯ , c n and how it results from drug effects on replication and death rates for different modes of action.

| Effect of monotherapy for different modes of action
We describe the net per-capita growth rate in the presence of treatment with a single drug (n = 1) by a sigmoidal Hill function using the formulation by Regoes et al. (2004): This dose-response relationship describes the pharmacodynamics of the drug.The susceptibility parameter zMIC W is the concentration at which the net growth rate is zero.In the limit of very high concentrations, the net growth rate converges to min (lim c→∞ (c) = min with  min < 0h −1 ).The Hill coefficient describes the steepness of the curve.
We assume that the Hill coefficient and the parameter min are the same for all cell types, which can be observed in experimentally measured dose-response curves (Chevereau et al., 2015;Das et al., 2020), but might not to be universally true.The cost of resistance affects the intrinsic replication rate, as discussed earlier, and resistance provides a benefit of increasing the wild-type parameter zMIC W by a factor res .The net growth rate of a resistant mutant is thus calculated by The drug effect on the net growth rate can be due to an increase in the death rate, a decrease in the replication rate or both.We refer to this as the 'mode of action' of a drug.Note that according to our definition, 'mode of action' thus does not relate to the specific mechanism of action or target site.For simplicity, we will drop the index indicating the cell type when describing the modes of action.

| Bactericidal drug that increases the kill rate independent of replication (CK)
The first mode of action that we consider are bactericidal drugs (indicated by 'CK') that increase the kill rate independent of replication and have no effect on the replication rate.An example for this drug type are aminoglycoside antibiotics, which are bactericidal and able to kill non-replicating cells (McCall et al., 2019;Mutschler et al., 2005).We thus have CK (c) = 0 and CK (c) = 0 + (c), where the function (c) describes the drug-induced killing.Comparing the net growth rate to Equation (1), we see that (c) = E(c) with E(c) given in Equations ( 2) or (3).We thus have CK (c) = 0 + E(c).

| Bacteriostatic drug that reduces the replication rate (S)
At the other end of the spectrum, we consider bacteriostatic drugs, such as chloramphenicol, tetracycline or macrolide antbiotics (Mutschler et al., 2005), that only affect the replication but not the death rate (indicated by 'S').The replication rate is reduced by a fac- The net growth rate is given by from which we can read off that For this mode of action, the net growth rate at high concentrations, min , can be at most − 0 , as the function (c) cannot be larger than one. (1) (2) . (3)

| Bactericidal drug that acts on replicating cells (CR)
Many bactericidal drugs only act on metabolically active cells, such as -lactam antibiotics, as well as drugs from other classes, for example ciprofloxacin or rifampicin (McCall et al., 2019).This is not captured by our CK drug, where killing is independent of the replication rate 0 (which is correlated with the metabolic activity).In particular, some bactericidal drugs act during replication itself, and we include this as a third mode of action.We describe the fraction of cell replications which result in cell death by a function (c) (with (c) ∈ 0, 1 ∀ c ).
The replication rate is thus reduced to CR (c) = 0 ⋅ (1 − (c)) and the death rate increased to CR (c) = 0 + (c) ⋅ 0 , which gives the net growth rate The function (c) is thus given by (c) = E(c) 2 and CR (c) = 0 + E(c) 2 .The lower limit of the parameter min is − 0 − 0 , as (c) ≤ 1. Table 1 provides a summary of the net growth rates and replication and death rates for all three modes of action.

| Modelling the combined effect of multiple antibiotics
Modelling the combined effect of multiple drugs requires choosing a null model for the effect of the drug combination in the absence of drug-drug interactions.The most common reference models are Bliss independence (Bliss, 1939) and Loewe additivity (Loewe, 1928;Loewe & Muischnek, 1926).Although the two approaches are often presented as 'rival approaches' (Greco et al., 1995), a glance at the original literature reveals that they are defined for different (complementary) situations-Bliss independence for drugs acting on different target sites and Loewe additivity for drugs with the same target.
In the following, we will give a brief overview of the history of the two concepts, before we describe how we implement Bliss independence and Loewe additivity in our model.

| Background on the history of Bliss independence and Loewe additivity
Bliss independence and Loewe additivity were developed as reference models for 'no interaction' to determine whether a combination of drugs can be classified as synergistic or antagonistic or, in other words, whether the combined effect of the drugs would be larger or smaller than expected under the assumption of zero interaction.
However, answering how exactly the combined effect in the absence of interactions can be calculated is not trivial.Let us have a look at the original definitions of the additivity models.
Before discussing the seminal publications by Loewe and Bliss (Bliss, 1939;Loewe, 1928;Loewe & Muischnek, 1926), we want to highlight the work by Wilhelm Frei, published more than 10 years before Loewe's first paper on drug additivity (and, just as Loewe's first paper, written in German).In 1913, Wilhelm Frei discusses, on a theoretical basis, how the joint effect of a drug combination could be predicted in order to compare the effect of combinations with monotherapy for a fixed total dose, varying the ratios of the drugs in combination (Frei, 1913): Let us assume we use two drugs A and B, which we apply in monotherapy at concentrations a and b , respectively.The concentrations are chosen in a way that the effects are the same (i.e.Frei (1913) for the visualisation of this example).He calls the additivity of doses 'Iso-Additivität' (iso-additivity) and the additivity of effects 'Hetero-Additivität' (hetero-additivity), which would be the same, for a linear relationship of the dose-response curves.He comments that iso-additivity is only possible for drugs with the same effect curves, but he does not mention any assumptions specifically for the target sites of the drugs.
Siegfried Walter Loewe mentioned the concept of additivity first in his work from 1926 (Loewe & Muischnek, 1926), explaining that two drugs would have what he calls a not-varied joint action if they behaved as if two doses of the same drug were combined, which would be 'additive'.He further explains that this can only be observed for two drugs with the exact same site of action.He uses an isobole diagram to explain the concept of additivity (see figure 1 in Loewe and Muischnek (1926)-additivity is given by line number II).An isobole diagram can be described as follows: consider a graph in which the x-axis displays the concentration of drug A and the yaxis of drug B. Both drugs lead in monotherapy to a certain effect e (6) TA B L E 1 Net growth rate, replication and death rate for the monotherapy of drugs with different modes of action.

Mode
).An isobole is a line which connects the points a and b, displaying all the points (all pairs of concentrations ã, b ) that result in combination in the effect e. Two drugs are additive if this line is linear, hence if it includes all points ã, b for which holds Equation ( 7) defines the isobole equation, which Loewe states 2 years later in Loewe (1928).Basically, under additivity, a portion of one drug can be replaced by a portion of another drug, without changing the effect, which is essentially 'iso-additivity' as described in Frei (1913).Having gotten aware of Frei's work, Loewe discusses both iso-and hetero-additivity in his work from 1928 and acknowledges that additivity can have multiple meanings: hetero-additivity would not result in a linear isobole and might in fact, look like an antagonism or synergism, which needs to be taken into account when investigating drug-drug interactions for specific combinations.
Chester Ittner Bliss defines 'independence' (he calls it 'independent joint action') through stochastic independence, specifically referring to drugs with different toxic actions (Bliss, 1939).The proportion of cells surviving the treatment with two drugs (drugs A and B) is given by the product of the proportions surviving the treatment with either one of the drugs: where p A and p B are the proportions of cells killed by the treatment with either drug A or drug B, respectively, and p A,B is the proportion not surviving the combined treatment.Hence, the two drugs act independently of each other.He additionally defines the concept of 'similar joint action', which is essentially (iso-)additivity, as described by Frei and Loewe, for two drugs targeting the same set of receptors.As Bliss independence is defined on the surviving fraction of a population, not on the net growth rate of cells, recent modelling studies often refer to the definition of Bliss independence on rates as given by Baeder et al. (2016).In this definition, the joint kill rate of the combination is given by the sum of the individual kill rates of the drugs, which results in an effect-addition, similar to the description by Frei under the name hetero-additivity.
Both Bliss and Loewe make a clear assumption on the target sites of the drugs for which their models apply.Nevertheless, the models are often seen as rival approaches in which one should universally be favoured over the other.Greco et al. (1995) et al. (2016).The authors simulated the combined drug effect of two drugs with a multi-hit model, explicitly simulating the drug-target binding.Bliss independence (defined on rates) was able to predict the simulated drug effect for drugs with different target sites and Loewe additivity for drugs with the exact same target site.

| The implementation of Bliss independence and Loewe additivity
Following the assumption for the sites of action of the drugs, Bliss independence can be implemented for all combinations of drugs differing in their dose-response curves or modes of action.For the joint effect of the drugs on the kill rate, we follow Baeder et al. (2016) who showed that Bliss independence translates into additivity of kill rates.For the effect on the replication rate, we assume stochastic independence, that is, probabilities for successful replication multiply.
Table 2 shows the resulting replication and death rates for combining two drugs with the same or different modes of action.We would like to point out that neither of these implementations leads to stochastic independence of survival probabilities (based on our derivation of extinction probabilities in Appendix C) and thus does not actually match the original definition by Bliss (1939); see Appendix B.
For Loewe additivity, as drugs act on the same target site, we consider only combinations of drugs with the same mode of action.
To calculate the joint drug effect, we follow the implementation of Lederer et al. (2018).For two drugs which have in monotherapy the net growth rate functions X,1 c 1 and X,2 c 2 for a type X, the net growth rate under the combination treatment can be calculated by Equations ( 9a) and (9b) lead to the same results only under the condition that the dose-response curves are the same (up to potentially different minimal inhibitory concentrations), which we assume in the following when considering Loewe additivity.
To compare different treatments (monotherapies and combinations), we scale the individual drug concentrations with the respective MICs and effectively keep the total drug concentration constant.More precisely, we choose the individual concentrations is the same for all treatments, irrespective of the number and types of drugs in the treatment.Drug ratios refer to the scaled concentrations.For example, for a two-drug combination with a 1:1 ratio, we have c 1 ∕ zMIC W,1 = c 2 ∕ zMIC W,2 .For notational simplicity and concreteness, we set zMIC W,i = zMIC W = 1 mg mL for all drugs in the following, but results apply more generally.
Under Loewe additivity, the net growth rate is independent of the number of drugs (displayed for CK drugs in Figure 1a,c).For Bliss independence (panels b and d), by contrast, the number of drugs matters.Two observations can be made.First, the total concentration at which the net growth rate is zero (the MIC of the combination) depends on the number of drugs and the maximum strength of the drug effect min (Figure 1a,c and Figures S1, S2).Interestingly, the increase in the concentration required for growth inhibition, which we see for drugs with = 2 and low min (i.e.| min | large), can be observed in experiments (Russ & Kishony, 2018).Unfortunately, the study did not consider drugs with shallow dose-response curves.
Therefore, we do not know if the decrease in the concentration needed for growth inhibition that we observe for those curves can be observed empirically as well.
The second observation from Figure 1 is that the growth rate in the limit of high concentrations decreases with an increasing number of drugs as we are summing up the kill rates of the drugs.It seems rather unlikely that the effect of multiple drugs is so much stronger than that of a single drug in isolation and that this effect increases indefinitely with the number of drugs, which shows that there are puzzles with translating Bliss independence to rates.
We come back to this in the Section 4. For combinations of drugs with the other modes of action, the maximum effect can also increase with the number of drugs, but the minimum growth rate of susceptible bacteria is limited to − 0 for bacteriostatic drugs and − 0 − 0 for bactericidal drugs that act during replication (see Figure S1).We can also observe a higher maximum effect for combinations than for the respective monotherapies if the drugs have different modes of action, especially also-depending on the parameters-for combinations of bacteriostatic and bactericidal CR drugs (see Figure S3).
Examples of dose-response curves (net growth rates) for the resistant cell types and for different drug combinations are shown in Figures S1-S8.

| The probability of treatment success
For the comparison of different treatments, we calculate the probability of treatment success, which we define as the extinction probability of a bacterial population of initially N 0 cells undergoing treatment.
TA B L E 2 Replication and death rates of a cell with type X for the combinations of two drugs with different modes of action assuming Bliss independence.
The functions X,i , X,i and X,i with i ∈ 1, 2 are derived from the action of drug i in monotherapy as given in Table 1.The parameters describing the drug effect are specific to each drug (i.e.i , min,i ).

F I G U R E 1
The net growth rate of susceptible bacteria for monotherapy and combinations with up to four drugs.The functions are displayed for the different additivity models (columns) and different Hill coefficients (rows).For Loewe additivity, all lines coincide.The parameter values are given in Table 3, except min was set to − 1h −1 .

| Exact calculation of the probability of treatment success
In our model, cells replicate and die independently from each other (see Section 2.1).Hence, each cell in the initial population founds a lineage with its own independent fate-extinction or growth.To calculate the extinction probability of the entire population, we need to multiply the extinction probabilities of all these independent lineages.For the treatment with n drugs, the probability of treatment success is thus given by where q X c 1 , c 2 , ⋯ , c n is the extinction probability of a cell lineage starting with a cell of type X (including the possibility of de novo mutations) and N 0,X cells of type X are present in the initial population.Importantly, even cells with a positive net growth rate might fail to establish a long-term lineage due to stochasticity at low cell numbers.The extinction probabilities q X are derived in Appendix C.

| Approximation of the probability of treatment success
Equation ( 10) allows us to compare the efficiencies of different treatments but we cannot see from it why a certain treatment is better than another.We therefore derive an approximation with which we can disentangle where the benefit of a treatment comes from, using an approach that is common in models of evolutionary rescue (see e.g.Bell & Collins, 2008;Czuppon et al., 2021;Martin et al., 2013;Orr & Unckless, 2008;Uecker et al., 2014).In the main text, we only give the derivation for monotherapy; the extension to combination treatment with two drugs is given in Appendix D.
We first focus on the de novo evolution of resistance.Two factors determine whether the population adapts or goes extinct: the number of resistant mutants appearing by mutation during treatment and the probability that these mutants establish, that is, do not suffer stochastic loss.New mutants appear at rate u where W(t) is the wild-type population size at time t.Given the large initial population size, we approximate the decline of the wild type deterministically: As the appearance of mutants follows a Poisson process, the number of mutants over the time course of the whole treatment is given by a Poisson distribution with mean Each of them establishes a long-term lineage with probability with q M being the extinction probability as derived in Appendix C.
The probability of treatment success, that is, the probability that no successful mutant appears over the time course of the treatment, is then given by the zeroth term of a Poisson distribution with mean N (tot) M p M est : We hence see that the probability of treatment success is determined by the product of the total number of new mutations appearing and their establishment probability.Single mutants in the standing genetic variation are eradicated with probability M and p M est are small, where the importance of minimising the establishment probability increases with the number of pre-existing mutants.By calculating these terms and comparing them among treatments, we can gain insights into whether a treatment is particularly successful because it prevents the appearance of new mutants or minimizes the establishment probability or both.

| RE SULTS
To identify the optimal conditions for antibiotic combinations, we will compare in this section the success probability of different combination treatments, including a comparison to monotherapy at the same total dose.We evaluate the treatment efficiency for several mutation probabilities towards resistance, ranging from low to high, and in the absence and presence of pre-existing mutants.We do so across a large range of concentrations.Absence of pre-existing mutations is unlikely, unless mutation probabilities are very low.Yet, it is instructive to consider this scenario, and it can be interpreted as preventing the de novo evolution of resistance.The parameter values used for the numerical analysis are displayed in Table 3.
Examples of drugs and drug combinations that are used to treat infections, together with key properties relevant to out study, are listed in Table S1.Whenever possible, we derive results for both additivity models (Bliss and Loewe) to evaluate if they lead to the same conclusion.
We first explore the role of the mode of action on the evolution of resistance in monotherapy.We then compare the efficiency of monotherapies with combinations and subsequently focus on combinations of drugs with different modes of action.Lastly, we extend the comparison to treatments at which the drugs are administered at different ratios and consider combinations with more than two drugs.

| Comparing monotherapies of drugs with different modes of action
Before considering multiple drugs in combination, we first compare drugs with different modes of action in monotherapy. (10) 3.1.1| Bacteriostatic drugs limit the de novo evolution of resistance better than bactericidal drugs if both have the same dose-response curve reveals that the bacteriostatic drug leads to a higher probability of treatment success than the bactericidal drugs and that the bactericidal drug acting during replication (CR) is better than a bactericidal drug which kills independent of replication (CK).To explain the difference in treatment efficiency, we compare the effect of the treatments on the number of mutants arising de novo and their establishment probability.Figure 2c shows that fewer mutants are generated during the treatment with the bacteriostatic drug compared to the bactericidal drugs (S < CR < CK).However, bactericidal drugs are better at preventing mutants from establishing, with the order being reversed, that is, CK < CR < S. A mathematical proof shows that these observations generalize beyond the numerical example (Appendix E): The more drug effect is allocated on reducing the replication rate versus increasing the death rate, the lower the number of new mutants, the higher the establishment probability and the higher the total probability of treatment success; that is, the effect on the number of mutants outweighs the effect on the establishment probability.
3.1.2| Bactericidal drugs are better than bacteriostatic drugs at preventing mutants from establishing, which can turn them into the better option when mutants pre-exist prior to treatment When mutants pre-exist in the population before the treatment, the effect on the establishment probability gains in importance.Figure 2b (right) shows that for drug concentrations above ca.2.5 ⋅ zMIC W , the bactericidal drug acting independent of replication leads to the highest probability of treatment success and the bacteriostatic drug to the lowest.At such high concentrations, de novo mutants are rare anyway, making it more important to prevent pre-existing mutants from establishing than to further reduce the number of new mutations.Note that establishment probabilities are generally high for the chosen parameter values such that even a single pre-existing mutant has a strong visible effect.For low concentrations, the number of new mutants is high (panel c), and efficiently decreasing the number of mutants by using a bacteriostatic drug is still the better strategy.
3.1.3| Bactericidal drugs can become better than bacteriostatic drugs for high enough bactericidal effects even in preventing the de novo evolution of resistance In reality, the value of min of the bactericidal drugs is most likely lower than that of the bacteriostatic drug; this especially holds for the Hill coefficient 0.5,1 or 2 Note: When multiple values are given, these were varied during the analysis.The pharmacodynamic parameter values were chosen in accordance to values measured from experimental data (Regoes et al., 2004).Values for the cost and benefits of mutations were reviewed by Igler et al. (2021); the values chosen here specifically rely on the work by Melnyk et al. (2015) and Spohn et al. (2019).
We consider a wide range of mutation probabilities, including values on the extreme high and extreme low end.Examples of these values can be found in Rodríguez-Rojas et al. ( 2014), Imhof and Schlötterer (2001), Oliver et al. (2004), McGrath et al. (2014) and Wang et al. (2001).Please note that we use the same zMIC W for all drugs.However, as the concentration range is always displayed in multiples of the zMIC W , the results are independent of the actual values, which could differ between drugs.

28
TA B L E 3 Parameter descriptions and their values used to generate the results.
CK drug.We therefore include such a drug comparison (Figure 2d).
In this comparison, the bacteriostatic drug becomes the worst option as the other two drugs lead to fewer new mutants (panels f).
At first sight, surprisingly, the mutant establishment probability is higher now for the bactericidal drugs than it was for the high value of min (i.e.| min | small).This can be understood by considering the dose-response curves: below the MIC, bactericidal drugs with a lower min have a higher net growth rate (panels a and d show this for the wild type, but the same applies to the mutant).
Overall, using a drug with a mode of action that strongly decreases the number of mutants, even at the cost of allowing resistant mutants to establish more easily than with other drugs, seems to be a good treatment strategy when resistance evolves de novo.When mutants are present in the initial population, using a drug that reduces the mutants' establishment probability most efficiently is better.

| Comparing monotherapies with two-drug combinations at equal drug ratios
We compare in this section the treatment efficiency of monotherapy with that of combination treatment with two drugs, given a constant total dose.For combination therapy, the total drug dose is thereby split equally between the two drugs (for unequal MICs, the drug doses are scaled with the respective zMIC W , see model description).
We will, for now, consider combinations in which all drugs of the comparison have the same dose-response curves and focus on the influence of the Hill coefficient .We mainly present results for bac-

| Combination therapy is always better than monotherapy for drugs with shallow doseresponse curves
Figure 3a-c shows that the probability of treatment success is always higher for combination treatment than for monotherapy when the dose-response curves are shallow (i.e.i ≤ 1) (the results for i = 1 are shown in Figure S9).This result is independent of the mutation probability and the pre-existence of mutants prior to treatment (Figure S10a-f).It holds both for combinations displaying Bliss independence and combinations displaying Loewe additivity.
Combinations with Bliss independence lead to a higher (or equal) probability of treatment success than drugs with Loewe additivity.
These observations are the same for all modes of action of the drugs (see Figures S11a-f and S12a-f).
What is causing the advantage of combination treatment over monotherapy?Dissecting the problem following the rationale of approximation ( 14) shows that for the chosen parameters, the number of new single mutants (for combination therapy, the sum of both single mutants) that appear during treatment are similar for both treatments, yet their establishment probabilities are much lower under combination therapy (see Figure S13).The establishment probability of the double mutants is similar to that of the single mutant in monotherapy but they are rarely generated.One can also understand this by considering the dose-response curves of all types (Figure S4): both susceptible cells and single mutants have lower net growth rates under combination therapy than under monotherapy.
For ≤ 1 , double mutants have similar or lower growth rates under combination therapy than under monotherapy (but as said before, they are rare).Combination therapy is thus better than monotherapy by decreasing the establishment probability of the single mutants and increasing the genetic barrier to resistance.The lower establishment probability of the single mutants is also intuitively clear beyond the numerical example: in monotherapy, they effectively experience a concentration that is lower by a factor (which we set to 28, see Table 3) than the one experienced by the wild type.In combination therapy, they still experience one half of the second drug to which they are susceptible.

| Monotherapy can sometimes be better than two-drug treatments for drugs with steep doseresponse curves
For a combination of drugs with steep dose-response curves, we can sometimes observe small concentration ranges in which monotherapy F I G U R E 3 Comparison of the probabilities of treatment success under monotherapy (n = 1) and a two-drug combination (n = 2).The results of the two-drug combination are shown for both additivity models, Bliss independence (pink) and Loewe additivity (green).Panels a-f display the results for bactericidal drugs acting independent of replication.Each column displays the comparison for different mutation probabilities (see panel header of the first row) and each row for different Hill coefficients (see values of in the header).Panels g-i display additional results for drugs with i = 2 and u i = 10 −11 (i.e. the parameter settings of panel d).Panel g considers the pre-existence of one single mutant per type in the initial population and panel h and i display the results for the other two modes of action.
leads to a higher probability of treatment success than combination therapy (Figure 3d-i).First, monotherapy can be slightly better than drug combinations displaying Bliss independence for concentrations slightly above the MIC (zMIC W ) of the wild type (panels d, g and h).
This can be explained by the increased MIC for combinations that we already saw in Figure 1 (see also Figure S4)-slightly above zMIC W , monotherapy can eradicate the wild type while combination therapy cannot.As briefly indicated above, whether this shift in MIC occurs or not, depends on the min of the drugs.For the parameter values used here, the shift in MIC occurs-opening up a small window in which monotherapy is superior-for both bactericidal drugs, but not for the bacteriostatic drugs (panels d, h and i; see Figures S4-S6).
Second, under Bliss independence, monotherapy leads to higher treatment success than combination therapy when drug concentrations are around res ⋅ zMIC W (the MIC of single mutants in monotherapy) and mutation probabilities are high.In this regime, concentrations in monotherapy are so high that even resistant types are controlled well.In contrast to drugs with shallow dose-response curves, double resistant types in combination therapy have a higher growth rate than single mutants in monotherapy at those concentrations and can still grow (see Figure S4 for the growth rates of all types).At the same time, they arise at a substantial frequency due to the high mutation probability.This effect is either absent or negligibly small for the other modes of action (Figures S11, S12).
When mutants pre-exist, monotherapy is usually much less effective, which further reduces the drug ranges in which monotherapy is better than combination treatment (see Figure 3g and Figure S10 for more results).An important exception to this are scenarios where mutation probabilities are high (which also implies many pre-existing mutants).In this case, monotherapy is superior to two-drug combinations with Bliss independence because types that are fully resistant to the respective treatment are better controlled by monotherapy than by combination therapy (Figure S10).
To conclude, we see that although there are some noteworthy exceptions, the two-drug combination seems to be overall a better treatment choice than monotherapy across different mutation probabilities and for all three modes of action, despite the decrease in the individual drug doses.

| Comparison of combinations of two drugs with different modes of action (50:50 ratio)
In this section, we discuss combinations of drugs which differ in their modes of action.For the comparison, we will again match the dose-response curves of the individual drugs as in Section 3.1, although-as pointed out abovemin is likely lower for bactericidal than for bacteriostatic drugs.The dose-response curves of the combination (i.e. the function c 1 , c 2 ), however, differ between the modes of action of the drugs (Figure S3).As in the previous section, the total drug concentration is split equally between the two drugs.We will first compare combinations of drugs with different modes of action to combinations of drugs with the same mode of action to evaluate if mixing drugs of different modes of action is generally beneficial.Then, we compare various combinations of drugs with different modes of action to identify which one is optimal.We show results for drugs with a Hill coefficient i = 1.Additional results for other parameter choices can be found in Appendix S1: Section S4.Varying the intrinsic replication and death rate and the Hill coefficients does not change the outcome of the comparisons.

| Combining drugs with different modes of action can sometimes be the best option by reducing both the number of mutants and the establishment probability
We find that using drugs with different modes of action can lead to a higher probability of treatment success than treatments using two drugs with the same mode of action.This is shown in Figure 4a-c: combining a bacteriostatic drug with a bactericidal drug acting independent of replication (S + CK) is better than using two bacteriostatic drugs (S + S) and at least as good as combining two CK drugs.
To see what makes this combination so successful, we can look at the number of mutants generated and their establishment probabilities; this is shown in panels d-f for a high mutation probability where the CK + S combination is better than both of the other combinations.
Combining two bactericidal drugs leads to the lowest establishment probability but the highest number of mutants, and combining two bacteriostatic drugs leads to the lowest number of mutants but the highest establishment probability (panels d and e).By suppressing both factors fairly well, the mixed treatment becomes the best treatment option.Hence, in this particular case, combining different modes of action is advantageous.However, a similar comparison with a bacteriostatic (S) and a bactericidal drug acting during replication (CR) shows that the mixed treatment (S + CR) never results in the highest probability of treatment success (Figure S14).

| Combining a bacteriostatic drug with a bactericidal drug can be both the best or worst design of a combination treatment
Figure 5 compares combinations of two drugs with different modes of action.The best treatment option, in the absence of pre-existent mutants, is for all mutation probabilities the combination of the bactericidal drug acting independent of replication (CK) with a bacteriostatic drug (S) (Figure 5a-c).When mutants pre-exist in the initial population, combining drugs with the two bactericidal modes of action (CK + CR) is sometimes slightly better (panels d-f).Combining the bactericidal drug acting during replication (CR) and the bacteriostatic drug (S) is always the least good option.Results are similar for drugs with Hill coefficients i = 0.5 and i = 2 with the noteworthy exception that the S + CR combination can be better than the CR + CK combination if i = 2 and the mutation probability high (Figures S15 and S16).All combination treatments, even the worst one, are better than the best monotherapy (see grey line in the plots).
We can again look at the two factors that determine resistance evolution-appearance and establishment of mutants (panels g-i).
Combining the bacteriostatic drug with the bactericidal drug acting independent of replication (CK + S) leads to the lowest number of mutants.This combination suppresses the establishment probabilities as much or almost as much as the combination of two bactericidal drugs, which has the lowest establishment probabilities.
For the combination of the bacteriostatic drug with the bactericidal drug acting during replication (S + CR) is bad in terms of establishment probabilities and-somewhat surprisingly-also in terms of number of mutants generated.The result does not change much when the first drug (the drug that is first written in the legend) has a higher maximum effect, that is, lower min (Figure S17).

F I G U R E 5
Comparison of the probabilities of treatment success for various combinations of drugs with different modes of action and for the best monotherapy.The first and second row show the probabilities of treatment success for three different mutation probabilities in the absence (a-c) and in the presence (d-f) of preexisting mutations.As reference, the best monotherapy is displayed in grey.The third row shows the components of the approximations of the probability of treatment success (Equations ( 14), (A16) and (A21)) for the highest mutation probability (u i = 10 −5 ).All drugs have the same dose-response curve with Hill coefficient is i = 1 and Comparison of the probabilities of treatment success for combinations of drugs with either the same (S + S and CK + CK) or different modes of action (S + CK).The upper row shows the probabilities of treatment success for three different mutation probabilities.The lower row shows the components of the approximations of the probability of treatment success (Equations ( 14), ( A16) and ( A21)) for the highest mutation probability (u i = 10 −5 ).
All drugs have the same dose-response curve with Hill coefficient is i = 1 and Overall, combining drugs with different modes of action can contain both the establishment probability of mutants and the number of mutants, resulting in a more efficient treatment than a combination that mainly influences one of the two factors.However, not all combinations of modes of action are equally good.Combining a bacteriostatic drug with a bactericidal drug killing independently of replication (S + CK) is the best (or very close to the best) treatment option across all comparisons.Combining a bacteriostatic drug with a bactericidal drug that acts during replication (S + CR) is one of the worst combinations.However, it is still much better than the best monotherapy.

| Comparison of two-drug combinations at unequal drug ratios
In the previous sections, we discussed comparisons of treatments in which drugs were either administered in monotherapy or in a two-drug combination in which the total dose was equally split between the two drugs.In this section, we want to extend the comparison of the previous sections to combinations in which the drugs are not administered at a 50:50 ratio.We will provide an overview of different examples to show how the drug characteristics can influence the outcome of combinations administered at uneven drug ratios.The results in this section are displayed for drugs with Bliss independence; more results can be found in Appendix S1: Section S5.

| The ideal drug ratio might vary across the concentration gradient and across mutation probabilities even when combining drugs with the same drug characteristics
For the combination of drugs with the same shallow dose-response curves ( i ≤ 1), we can see, independent of the total drug dose, the mode of action, the mutation probability and the preexistence of resistance mutations, that the two drugs should be given at an equal ratio to maximize the probability of treatment success (see panels a-c in Figure 6  We can further ask which drug ratio minimizes the concentration needed to have an almost certain treatment success, which we define here as P success ≥ 99 % (the white area in the panels of Figure 6).For shallow dose-response curves, the minimum concentration is lowest for a 50:50 ratio (the boundary between the coloured and white areas is concave).For steep dose-response curves, a 50:50 is optimal for low mutation probabilities.In contrast, for high mutation probabilities, monotherapy minimizes the concentration required for almost certain treatment success, and in combination treatment, much higher concentrations would be needed to achieve a treatment success of 99 %.Yet, at the concentration where monotherapy gives a 99 % probability of treatment success, administering drugs at a 50:50 ratio still leads to When mutants pre-exist, choosing the drug ratio that minimizes the concentration required for an almost certain treatment success (the boundary between the coloured and white areas) becomes more important as the success probability changes more drastically with the drug ratio (Figure S20).
For drugs combining according to Loewe additivity, a 50:50 ratio always minimizes the concentration required for almost certain treatment success (Figure S21).

| The optimal drug ratio varies when combining drugs with different drug characteristics
When combining drugs with different Hill coefficients (Figure 6g-i) or different modes of action (panels j-l), it can be beneficial to administer two drugs at unequal doses; this is especially true if mutation probabilities are high.Unlike for a combination of drugs with the same characteristics (as considered in the previous section), the ratio that minimizes the concentration required for almost certain treatment success is not either monotherapy or 50:50, but can be something in between.Yet, in many cases, the optimal ratio both in terms of maximising P success for a given concentration and in terms of minimising the concentration required for P success ≥ 99 % is still 50:50 or close to it.In cases where it is not, the mistake made by administering drugs at equal doses seems to be small (although it is unclear what should be considered 'small' here, given that the lives of patients may depend on it).In contrast, in cases where an unequal ratio is optimal (e.g.panel l), deviating from the optimal ratio in the other direction (i.e.giving even more of the high-dosed drug) can be detrimental, that is, strongly reducing the probability of treatment success.
Overall, our results show that an equal distribution of the drug concentration among two drugs in combination might not always be the ideal drug ratio.While it is mostly the best strategy when the two drugs have the same drug characteristics and generally for low mutation probabilities, an unequal ratio is often beneficial when the two drugs vary in their Hill coefficients or modes of action, and the mutation probability is high.Nevertheless, the examples considered here suggest that, when information is limited, a 50:50 ratio is a good choice as it is usually not much worse than the best strategy (an exception are drugs with steep dose response curves and high mutation probabilities, as discussed before).

| Comparison of combinations with more than two drugs
So far, we have focused on combinations with only two drugs, but what happens if we increase the number of drugs?In this section, we consider combinations of up to seven drugs to generalize our results from Section 3.2.All drugs in combination have the same dose-response curve and a bactericidal mode of action, killing independent of replication (CK).The total drug dose is split equally among them.We display the results for different Hill coefficients.
As more drugs are expected to be especially important when resistance pre-exists, we show those results in the main text.Varying the modes of action does not change the results discussed in the following (Figures S22 and S23).S24).Yet, from a certain number of drugs, the benefit of adding more drugs becomes small.
The same is found if we only account for the de novo evolution of resistance, except that a two-drug treatment is then still much better than monotherapy even at high mutation probabilities (Figure S25).
The increase in treatment efficiency when increasing the number of drugs can be explained similarly to the results comparing the two-drug combination and monotherapy (see Section 3.2).For each type, the net growth rate decreases with the number of drugs (see Figure S8 for the net growth rates of susceptible cells, single and F I G U R E 6 Probability of treatment success for combinations of two drugs administered at unequal ratios.The x -axis gives the total drug concentration, the yaxis the proportion of the first drug, and the colour indicates the probability of treatment success.The white area indicates almost certain treatment success (P success ≥ 99 %).The red marker indicates the ratios which maximize the probability of treatment success for a specific concentration.The optimal ratio is not shown for cases in which all treatments are equally bad, or once the first treatment has reached a success probability of at least 99 %.The figure provides an overview over specific cases (see row titles), more results can be found in Appendix S1: Section S5. Results are shown for different mutation probabilities, as displayed in the column header.Panels a-f show results for combination of drug with the same Hill coefficient (either i = 1 or i = 2).Panels g-i display results for a combination of drugs with different Hill coefficients ( 1 = 1, 2 = 0.5) and panels j-l for drugs with different modes of action (both with min = − 0.1h −1 , but the same Hill coefficient i = 2).
double mutants).The more drugs, the more mutations are required for resistance and the less likely resistant types are generated.
3.5.2| Bliss Independence and Loewe additivity disagree: more is not always better when doseresponse curves are steep While the results for drugs displaying Loewe additivity do not change for drugs with steep dose-response curves, Figure 7d-f show that for Bliss independence, more drugs are not always better.This can be explained in the same way as for comparing twodrug treatment and monotherapy in Section 3.2: The wild type gets better suppressed by treatments with fewer drugs.The same holds for the type resistant to the respective treatment (single mutants get better suppressed by monotherapy than double mutants by two-drug treatments, etc.).However, any resistant type is controlled better by adding another drug (e.g.single mutants are better controlled by multi-drug treatment than by monotherapy).
Which effect dominates varies with the drug concentration, the mutation probability during treatment, and the pre-existence of mutants.Accordingly, the optimal number of drugs depends on the conditions.As already observed in Section 3.2, two-drug treatment performs worse than monotherapy if the mutation probability is high and pre-existence is considered (panel f).Here, adding a third drug provides a huge benefit as-while double mutants are frequent-triple mutants are not.
As a rule of thumb we can formulate that increasing the number of drugs increases the probability of treatment success up to some point, but then the benefit decreases or-in some cases-even turns into a disadvantage.Overall, our results suggest that, when mutants appear frequently and two-drug treatments are inefficient, combining three drugs might already provide a sufficiently large benefit compared to combinations with two drugs.
F I G U R E 7 Probability of treatment success for combination with up to seven drugs for Bliss independence and Loewe additivity, accounting for preexisting resistance mutations.The results are displayed for different mutation probabilities (columns) and different Hill coefficients (rows).All drugs are bactericidal, killing independent of cell replication (CK).

| DISCUSS ION
How do we best design combination treatments to limit the evolution of drug resistance?To answer this question, we set up a stochastic pharmacodynamic model for treatment comparison and explored how drug characteristics and dosing strategies affect resistance evolution.To understand where the benefits of a certain treatment come from, we look at two factors influencing the evolution of resistance: the number of mutants generated de novo and the establishment probability of the mutants, that is, the treatments' ability to prevent mutants from arising and to control the growth of mutants once present.

| The effect of the mode of action of drugs on the evolution of resistance
For monotherapy, we found that bacteriostatic drugs are better at limiting the de novo evolution of resistance than bactericidal drugs, provided the effect on the net growth rate is the same.Under this condition, bacteriostatic drugs are better at preventing mutants from arising than bactericidal drugs but result in a higher establishment probability.Bactericidal drugs only become advantageous when mutants pre-exist or for large enough bactericidal effects.
Only recently, stochastic models started to compare how the type of drug effect-reducing the replication rate or increasing the death rate-affects the evolution of resistance (Czuppon et al., 2023;Marrec & Bitbol, 2020;Raatz & Traulsen, 2023).In a pharmacody-  2023), unlike us, include strain competition (assuming that non-replicating cells are still metabolically active); this means that growth of resistant strains is suppressed by competition with the sensitive strain for treatments with bacteriostatic drugs.Whether a drug affects the replication or the death rate is also of importance in more complicated ecological or genetic settings, for example, it affects the rate of resistance evolution in fluctuating environments (Marrec & Bitbol, 2020) or the number of lineages that evolve during an adaptive walk through the trait space (Raatz & Traulsen, 2023).
Overall, the results highlight the importance of differentiating between drug effects on birth and death rates and accounting for stochastic population dynamics rather than assuming a deterministic net growth rate in which bacteriostatic or bactericidal drugs are implemented by simply varying the size of the drug effect (i.e.min ).Coates et al. (2018) demonstrate that the extinction probability predicted by a birth-death process (with experimentally measured rates) is consistent with the experimentally determined probability of stochastic extinction of susceptible cells (basically 'proving' our Equation (A11) but for the wild type).This nicely shows that the mathematical approach of birth-death processes is suitable to describe stochastic bacterial dynamics.Extending the work by Coates et al. (2018) for more drugs and environmental conditions would help to improve model implementations and predictions on the efficiency of treatments using drugs with different modes of action.

| Comparing monotherapy with a two-drug combination
In most cases, two-drug combinations, administered at equal drug ratios, are better than monotherapy.Exceptions are found for drugs with steep dose-response curves (i.e. with a Hill coefficient larger than one) that combine according to Bliss independence.For those drugs, monotherapy controls susceptible and fully resistant types better than the combination and is therefore sometimes the better treatment choice, especially when mutation rates are high and mutants pre-exist.A pharmacokinetic-pharmacodynamic model for sequential therapy, where drug doses overlap at each administration, previously found that the same effect turns the maximum cycling frequency sub-optimal for such drugs (Nyhoegen & Uecker, 2023).
In vitro studies generally confirm that combination treatment can reduce resistance evolution compared to monotherapy (Angst et al., 2021;Barbosa et al., 2018;Hegreness et al., 2008;Jahn et al., 2021;Munck et al., 2014), but show that the ability to reduce resistance evolution is influenced by drug-drug interactions (Barbosa et al., 2018;Hegreness et al., 2008) or collateral effects (Barbosa et al., 2018;Munck et al., 2014).In clinical studies, resistance evolution is rarely considered as an endpoint.Among the few studies that do so, some observed a decrease in resistance evolution when using antibiotic combinations rather than monotherapies; others, however, find no difference, leaving the comparison inconclusive for the clinic (Tamma et al., 2012).Similarly, a recent metaanalysis of 29 studies assessing the effect of combination therapy on the within-patient evolution of resistance in comparison to monotherapy did not identify a general effect (Siedentop et al., 2023).
However, the authors further showed that most of the clinical trails did not possess sufficient statistical power, highlighting the need for much larger clinical studies to fill this important knowledge gap.
Mathematical models which incorporate important aspects of patient treatment such as the pharmacokinetics of the drugs might help to understand how our findings and those from in vitro studies translate to the clinic.
We further extended our comparison to combinations of drugs at unequal drug ratios.While equal drug ratios are most often beneficial for low mutation rates, unequal drug ratios can be optimal at high mutation rates, especially when combining drugs with different characteristics.However, a 50:50 ratio is usually not much worse.To our knowledge, the effect of different drug ratios on the evolution of resistance has not yet been systematically studied.We only considered a limited set of examples here, and an in-depth exploration of the parameter space is required to further assess the importance of the drug ratio.As it is impossible to achieve a specific drug ratio within the patient, such an analysis needs to consider the errors made by deviating from the optimal ratio.

| Combinations of drugs with different modes of action
While combination treatments, irrespective of the modes of action, are usually better at limiting resistance evolution than monotherapy, not all combinations are equally good.For all cases considered here, the combination of a bacteriostatic drug with a bactericidal drug acting independent of replication was the best or one of the best choices.Combining instead a bacteriostatic drug with a bactericidal drug acting during replication led in most cases to the lowest success probability among the combination treatments included in the comparison (but the combination was still better than monotherapy).
Whether bacteriostatic and bactericidal drugs should be combined, is also of interest simply with respect to their potency to clear susceptible cells.Coates et al. (2018) showed in vitro that adding a bacteriostatic drug to a bactericidal drug that acts on non-replicating cells can strongly increase stochastic extinction of susceptible cells.
For bactericidal drugs that are only effective against replicating cells, in contrast, adding a bacteriostatic usually had no or even an adverse effect.Based on Loewe additivity, many combinations of bacteriostatic and bactericidal drugs are furthermore classified as antagonistic (Chandrasekaran et al., 2016;Ocampo et al., 2014) or even suppressive (Lázár et al., 2022).However, Ocampo et al. (2014) also found synergistic combinations of bacteriostatic and bactericidal drugs, especially for combinations including aminoglycosides, which have been found to kill non-replicating cells (McCall et al., 2019).
Overall, combining bacteriostatic and bactericidal drugs can be a good strategy to limit the evolution of resistance.However, choosing such a combination requires a good understanding of how the drugs affect cellular processes and how these respective processes influence each other.

| Combinations with more than two drugs
Increasing the number of drugs often increases the treatment efficiency.However, the benefit of adding another drug decreases with an increasing number of drugs.For drugs displaying Bliss independence, an increase in the number of drugs can even decrease the treatment efficiency when the drugs have steep dose-response curves, as already discussed for two-drug treatment.In this case, our results suggest using a two-drug combination (or even monotherapy) for very low and a combination of three or four drugs (depending on the mode of action) for large mutation probabilities.Increasing the number of drugs is a beneficial strategy for rapidly evolving strains.The current guidelines for the treatment of M. tuberculosis suggest using a four-drug combination in the first 6 months of treatment (WHO, 2022).Such multi-drug combinations are common for other rapidly evolving pathogens such as HIV and Malaria (Goldberg et al., 2012).The use of more than two drugs can even be necessary when only a small fraction of the cells in the population have an elevated mutation rate (i.e.mutator-strains): in a recent study combining laboratory experiments and modelling, Gifford et al. (2023) observed resistance evolution against the combination of rifampicin and nalidixic acid for a population of E. coli in which only 6 % of the bacteria had an elevated mutation rate compared to the wild type, while multi-drug resistance did not evolve in the absence of the mutator strain.
To provide a baseline comparison, we decided not to include drug-drug interactions in our analysis.Nevertheless, increasing the number of drugs can increase the frequency of drug-drug interactions (Katzir et al., 2019;Tekin et al., 2018).Hence, accounting for interactions-including possibly higher order drug interactions-is a relevant next step.

| Contrasting the results obtained for Loewe additivity and Bliss independence
Whenever possible, we included results obtained for both additivity models, Loewe additivity and Bliss independence.Generally, the conclusions drawn from the two models are similar.However, for drugs with steep dose-response curves, the conclusions sometimes diverge due to the shift in MIC with Bliss independence, as discussed above.The need for higher drug doses with an increasing number of drugs aligns with findings from in vitro experiments (Russ & Kishony, 2018).This study concludes that Bliss independence can better predict the effect of multi-drug combinations than Loewe additivity (see also Katzir et al., 2019).Given the assumptions on the target sites made by the two models, this could be explained by the decreasing probability that all drugs target the same cellular site with an increasing number of drugs.
As mentioned in the Section 2, there are puzzles with implementing Bliss independence on rates.For bactericidal drugs, the common approach of summing up the death rates leads to an implausible high maximum effect, especially for large numbers of drugs.This approach as well as our implementations for the other modes of action further actually do not align with the original definition by Bliss (1939).How to translate Bliss independence to replication and death rates is not trivial and raises fundamental questions that we lack the space to address here.Modelling the drug effect explicitly through drug-target interactions (extend-

| Limitations and extensions
A caveat to our study is that-as we focused on a detailed understanding of the observed patterns-we only considered a limited set of parameters.While we varied several important parameters to test for robustness and many explanations seem intuitive, we cannot rule out that there may possibly be regions in the empirically relevant parameter range where results deviate.A study that focuses on a parameter sensitivity analysis would therefore be a highly valuable complement.
Apart from this general caveat, one of the major limitations of our work is that we only considered one resistance mutation for each drug, while in reality there are several or many mutations that can confer various levels of resistance.Generally, more mutations confer resistance to low than to high drug concentrations, which could impair the efficiency of multidrug low-dose treatments.Single resistance further evolved in our model with the same probability and resulted in the same increase in MIC for all drugs.The optimisation with respect to the drug ratio might become more critical if resistance to one drug is more frequent or stronger than to the other.We did not include cross-resistance, which-depending on its strength-might make combination treatment lose its benefit over monotherapy.
We only considered the acquisition of resistance through chromosomal mutations.Yet, resistance genes are often located on (Carattoli, 2013).Pathogenic bacteria could either carry a resistance plasmid to begin with or acquire it through horizontal gene transfer from commensal bacteria (Francino, 2016;San Millan, 2018) or other pathogens during a polymicrobial infection (Orazi & O'Toole, 2019).It is reasonable to assume that, similar to mutations, the acquisition rate is proportional to the population size; for example, conjugative transfer is often modeled by a massaction kinetics term (Levin et al., 1979).However, the rate depends on the dynamics of the donor population, which might be affected by the treatment too or interact with the recipient population.For cells whose growth is inhibited by a bacteriostatic drug, resistance on a newly acquired plasmid might not be expressed as for example assumed in the model by Willms et al. (2006).Importantly, the acquisition of a multidrug resistance plasmid would lead to resistance to multiple drugs in one step.There are numerous further important differences to chromosomal resistance, for example, the risk of segregative loss of incompatible plasmids (Cullum & Broda, 1979;Ishii et al., 1978;Novick & Hoppensteadt, 1978) or the effect of conjugation on the establishment of the resistance plasmid (Geoffroy & Uecker, 2023;Novozhilov et al., 2005;Tazzyman & Bonhoeffer, 2013, 2014).Our model is thus most relevant for treatment with antibiotics where resistance relies on mutations and cannot be acquired through plasmids, but it also applies to other antibiotics as long as no resistance plasmid is available to enter the population.For reviews on modelling plasmid dynamics and an introduction to plasmid biology for modellers with modeling examples, we refer to Hernández-Beltrán et al. ( 2021) and Dewan and Uecker (2023).
Concerning the modes of action, it has been criticized that drugs cannot be strictly classified as either bacteriostatic or bactericidal (Pankey & Sabath, 2004;Wald-Dickler et al., 2018) and that furthermore, the static or cidal activity might depend on the antibiotic concentration (Pankey & Sabath, 2004).In a model that accounts for that, the drug would thus need to independently act on both rates, possibly with primarily bacteriostatic effects at low and primarily bactericidal effects at high concentrations.We also neglect many aspects of bacterial growth, especially competition, which has repeatedly been found to affect the evolution of single and multidrug resistance (e.g.Berríos-Caro et al., 2021;Czuppon et al., 2023;Pena-Miller et al., 2013).
We restricted our analysis to constant drug concentrations.In the patient's body, the drug concentration declines over time due to pharmacokinetic processes (Levison & Levison, 2009), which could reduce the efficiency of a combination compared to monotherapy when administered at the same total drug dose.In addition, the pharmacokinetics of the drugs could vary, for example, in the elimination or penetration rate; hence, besides the concentration, the drug ratio could change with time.In extreme cases, this might lead to temporal or spatial monotherapy, which can facilitate resistance evolution (Moreno-Gamez et al., 2015).Last, it would be important to also consider clinically relevant endpoints other than resistance evolution, such as the time until the infection is cleared.Here could lay a potential advantage of bactericidal over bacteriostatic drugs, as death rates can be much higher and hence eradication faster.

| CON CLUS ION
With our study, we step by step developed a set of guiding principles of how to design combination treatments.We hope that it is the starting point for future studies that assess the robustness of the conclusions for multiple mutations with different resistance levels, unequal mutation probabilities for the two drugs, pharmacokinetic processes and other factors, and extend the model by drug-drug interactions.Our study further highlights the need to better link the traditional phenomenological definitions of drug independence at the population level and the microscopic dynamics at the level of cells.A model that is consistent between the two levels would generate a common framework for mathematical modelling, experiments on single-cell dynamics and traditional experimental measurements at the population level and would thus allow for a much better integration of modelling work and empirical results.
probability of a type can always be determined recursively from those with more mutations.
We now proceed to give the formal derivation of the extinction probabilities q X (c) founded by a single cell of type X.For notational simplicity we drop the dependency on the concentration c.To calculate q X , we consider the probability generating function corresponding to the 'offspring' distribution in the multi-type Galton-Watson process.For a treatment with n drugs, there are 2 n possible cell types in the population, and the probability generating function f(s) is thus a vector of length 2 n , and the argument s is a vector of the same length.
The component of f(s) that corresponds to cell type X is given by where s X and s Y are the components of s corresponding to types X and Y, u X→Y is the probability that a daughter cell of type X is of type Y and U = ∑ Y u X→Y.The probability u X→Y is zero unless type Y carries additional mutations with respect to type X.
As an example, we explicitly show the probability generating function for n = 2.In a treatment with two drugs, there are four possible types-the fully susceptible wild type (W), the fully resistant double mutant (M 1,2 ) and the two single resistant types (M 1 and M 2 ).
The probability generating function f(s) The vector of extinction probabilities q (e.g.q = q W , q M 1 , q 2 , q M 1,2 for treatment with two drugs) is the smallest fix point of the probability generation function (Sewastjanow, 1974, p. 115).One can convince oneself that the fix point equation q = f(q) corresponds to the equations that one obtains with the intuitive arguments in the main text.

A PPEN D I X D Approximation of the probability of treatment success for a twodrug combination
In this section, we derive the approximation of P success for a two-drug combination (drug 1 and drug 2), corresponding to approximation (14) for monotherapy.Each term-number of new mutants and their establishment probabilities-needs to be calculated for each type that is resistant at a given concentration.
We model the dynamics of all types that have a negative growth rate at that concentration deterministically and use this to obtain the number of new mutants resistant to the treatment.We then calculate their establishment probabilities without considering further mutations; that is, in concentration ranges where a single mutant has a positive net growth rate, we ignore the possibility that the single mutant might go extinct but generate a doubleresistant type before doing so.We only consider two-step rescue via single mutants when single mutants cannot grow themselves.
Hence, we split the concentration range into two parts.The first range includes all concentrations at which the single resistant types have a positive growth rate.For simplicity, we assume here that this range is the same for both single resistant types, which does not need to be true generally (in those cases the concentration range needs to be split into more parts).We formulate the establishment probabilities based on rates which follows from the transformation in Equation (A9).
For the first drug range, the mean numbers of mutants generated from the susceptible type (mean of a Poisson distribution) and their establishment probabilities are-nearly identical to the approximations for monotherapy-given by The probability of treatment success is then approximated by In the second drug range, we need to additionally consider the number of double mutants generated by the single resistant types .We are now interested in the p.g.f. of the random variable Z = ∑ Y j=1 X j , which is obtained by We further obtain for the probability that Z takes the value 0: where the approximation assumes that est is small (which is a very reasonable assumption).
We now still need to account for (1) the same pathway to resistance via single mutants resistant to drug 2 and (2) double mutants that get directly generated from wild-type cells by two mutations.
Setting we obtain for the probability of treatment success Given the low mutation probabilities, terms 1 − u 1 or 1 − 2 can be omitted as they are ca. 1.

A PPEN D I X E
Targeting the replication rate, the death rate or both: A comparison of monotherapies In Section 3.1 in the main text, we investigated the effect of the mode of action of drugs in monotherapy on the probability of treatment success.We observed that bacteriostatic drugs (S) are better than bactericidal drugs (CK and CR) and bactericidal drugs acting during replication (CR) are better than bactericidal drugs killing independent of replication (CK), given the dose-response curves are the same.In the following we will prove this mathematically.
Consider a drug where a fraction of the drug effect E(c) is allocated to decreasing the replication rate and a fraction (1 − ) of the drug effect is allocated to increasing the death rate, that is, = 0 − E(c) and = 0 + (1 − )E(c).For the S drug, we have = 1 ; for the CK drug, we have = 0; and for the CR drug, we have = 1 2 (see Section 2.2).
For such a drug, the total number of new mutants that get generated during population decline is given by (We are omitting the term (1 − u) here as it is small.)This is a decreasing function in , that is, the more the drug acts on the replication rate, the fewer mutations are generated.
For the establishment probability of mutants, we have which is an increasing function in ; that is, the more the drug acts on increasing the death rate, the lower the establishment probability of mutants.
What about the product of the two?
For the derivative of the function f( ), we obtain

It holds
In the regime that we are interested in-the wild type has a negative and the mutant a positive net growth rate-it holds As we assume that M 0 ≤ W 0 , it must especially hold that which shows that f � () < 0. This in turn means that the product M ⋅ p est decreases with .The smaller the product, the larger the probability of treatment success.Thus, the more drug effect is allocated to reducing the birth rate, the higher the probability of treatment success.This especially proves the order S > CR > CK.
E A (a) = E B (b)).All potential combinations of drugs A and B (with concentrations (p ⋅ a, (1 − p) ⋅ b) ∀ p ∈ (0, 1) ) could either result in the exact same effect as in monotherapy(i.e.E A (a) = E B (b) = E A,B (p ⋅ a, (1 − p) ⋅ b) )or in a larger (or smaller) effect.He illustrates this with the example of adding 1 6 b to the dose 5 6 a. Adding the second drug could, on the one hand, behave as if we would add just another fraction 1 E A (a) ), hence adding simply the doses of the two drugs together.On the other hand, adding 1 6 b could result in an effect addition with E Awhich can be different from E A (a), depending on the dose-response relationships of the drugs (see figure8 in

First
, we assume that the drugs have the same effect on the net growth rate, shown in Figure 2a.Panel b displays the probability of treatment success for a fully susceptible population and a population with one pre-existing mutant.Let us focus now on the evolution of resistance in the susceptible population.Panel b (left) Mutation probability, resistance to drug i 10 −11 , 10 −9 or 10 −5 zMIC W Wild-type MIC 1 μg/mL res Benefit of mutation in terms of resistance i Cost per mutation 0.1 min Growth rate at very high concentrations: tericidal drugs and only show specific cases for drugs with different modes of action as the results are otherwise the same.All results are shown for both models of drug additivity(Loewe and Bliss).We also varied the intrinsic replication and death rate and min .The results F I G U R E 2 Comparison of monotherapies with drugs with different modes of action.Panels a-c show the results for a comparison of drugs with the same dose-response curve ( min = − 0.065h −1 ), and panels d-f show the results for drugs with different doseresponse curves (differing in min ).Panels a and d display the respective doseresponse curves.Panels b compare the probabilities of treatment success for the three drugs in the absence and presence of pre-existing mutants.Panels e show the probabilities of treatment success for different mutation probabilities.Panels c and f show the terms used in the approximation of the probability of treatment success (Equations (13) and (12) and their product).The results are shown for = 1 and u = 10 −11 (except for panel e (right)).
, as the change in parameter values did not affect the outcome of the comparison.
and FiguresS18-S20, for drugs with the other two modes of action and in the presence of preexisting mutants).For drugs with steep dose-response curves, monotherapy is advantageous at low concentrations for low mutation probabilities (panels d and e) and high concentrations for high mutation probabilities (panel f).Both results-those for shallow and those for steep dose-response curves-align with our comparison of monotherapies versus two-drug combinations in Section 3.2.In addition to the explanations provided there, we show in FigureS7for all types the ranges where net growth rates are positive/negative.
3.5.1 | Bliss Independence and Loewe additivity agree: more is better when dose-response curves are shallow Considering drugs whose joint effect is described by Bliss independence leads to the same conclusion as for drugs with Loewe additivity if the drugs' Hill coefficients i ≤ 1: Increasing the number of drugs increases the probability of treatment success, even though the individual drug concentrations are reduced more with every drug added (Figure 7a-c and g-i and Figure namic model similar to ours, Czuppon et al. (2023) observe higher establishment probabilities of resistant mutants for treatments with bactericidal than with bacteriostatic drugs, which is opposite to what we find.The reason for this difference is that Czuppon et al. ( ing, for example, the work byBaeder et al. (2016)) might bring some clarity.Alternatively, modellers could work with specific experimentally measured dose response curves without caring how the combination is classified(Ankomah et al., 2013).For stochastic models, comparisons of time kill curves of drug combinations and monotherapies (such as inAnkomah et al. (2013) andYu et al. (2016)) would need to be extended by measurements of the replication and death rates of cells.

(
which are themselves doomed to extinction in this range).The number of single mutants generated by the wild type is again Poisson distributed by the mean given by the same expression as before (Equation (A15)).We further describe the decay of theses single mutants deterministically as well (although they are rare), and the number of double mutants generated from each of these single mutants is again Poisson distributed with means What is the probability distribution of the total number of double mutants generated in this way?Let us denote by Y the random variable describing the number of mutants resistant to drug 1 generated from the wild type.Y is Poisson distributed with mean N M 1 , and the p.g.f. is thus given by h Y (s) = e −N M ⋅(1−s) .Let us denote by X j the number of successful double mutants generated from the jth single mutant.X j is Poisson distributed with mean ÑM 1 →M 1.g.f. of is thus given by h X (s) = e ),M 1 est (c) = 1 − 0 + (1 − )E M (c) 0 − E M (c) >  W 0 −  0 − E W (c) Loewe additivity over Bliss independence.However, this brief excursion into the literature shows that neither of the models is universally applicable and rather complements the other than opposes it.Both models were defined for specific assumptions on the modes of action.This is further supported by the work of Baeder